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Abstract 

Lyapunov exponents (LEs) are key indicators of chaos in dynamical systems. In 
general relativity the classical definition of LE meets difficulty because it is not 
coordinate invariant and spacetime coordinates lose their physical meaning as in 
Newtonian dynamics. We propose a new definition of relativistic LE and give its 
algorithm in any coordinate system, which represents the observed changing law 
of the space separation between two neighboring particles (an "observer" and a 
"neighbor"), and is truly coordinate invariant in a curved spacetime. 
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Chaos is a popular phenomenon in dynamical systems. One of its main features 
is the exponential sensitivity on small variations of initial conditions. The 
exhibition of chaos in the motion of Pluto makes it particularly attractive for 
scientists to investigate the dynamical behavior of the solar systemfl]. 

In Newtonian mechanics, Lyapunov exponents (LEs), as a key index for mea- 
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suring chaos in a dynamical system, can be calculated numerically with ei- 
ther the variational method [2] or the two-particle method [3], and sometimes 
a mixture of the two [4]. The former is more rigorous but one has to derive the 
variational equations of the dynamical equations of the system and integrate 
them numerically with the dynamical equations together. The latter is less 
cumbersome especially when one wants to compute the maximum LE only, 
which is the main index for chaos. We will discuss the extension of the two 
particle method to relativistic models in this letter. 

Let q(t) and q[t) be the position and velocity vector of a dynamical system. 
As a set of initial conditions is selected randomly and the corresponding tra- 
jectories are restricted in a compact region in the phase space, the classical 
definition of the maximum LE is 



where the separation d between two neighboring trajectories is required to be 
sufficiently small so that the deviation vector (Aq, Aq) can be regarded as a 
good approximation of a tangent vector, and the distance d(t) at time t is of 
the form 



It is the Euclidian distance in the phase space, q and q are in different di- 
mension, so one has to carefully choose their units to assure that both terms 
in the expression of d(t) are important. We suggest computing the LE in the 
configuration space rather than in the phase space, that is, 
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We argue that both Ajv and A^y are effective in detecting chaotic behavior of 
orbits since Aq(£) = d/dt(Aq(t)) and they should have the same Lyapunov 
exponents. 

In general relativity the above definition and algorithm are questionable. First, 
there is no unified time for all the reference systems. Secondly, the separation 
of time and space of the 4-dimensional spacetime is different for different 
observers. Furthermore, time and space coordinates usually play book-keeping 
only for events and are not necessary with a physical meaning. Consequently, 
one would get different values of A^ and A^ in different coordinate systems. 

Up to now the references for exploring chaos in relativistic models have been 
mainly interested in studying the dynamical behavior of black hole systems [5- 
8] and the mixmaster cosmology [9-13]. Most of them adopt the classical def- 
inition of LE. Here we will abandon the variational method in the case of 
general relativity for it is rather cumbersome to derive variational equations. 
Actually there exist general expressions for geodesic deviation equations [5] 
but one has to derive the complicated curvature tensor. Furthermore, in many 
cases a particle does not follow a geodesic. For instance, a spinning particle in 
Schwarzschild spacetime doesn't move along a geodesic[5]. Therefore, we will 
concentrate on constructing a revised edition of the two-particle method. The 
classical algorithm of LE lacks invariance in general relativity, which depends 
on a coordinate gauge and even could bring spurious chaos in some coor- 
dinate systems. The chaotic behavior in the mixmaster cosmology[9-13] has 
been debated for decade or so. Using different time parameterization [14,15], 
one would get distinct values of A at. Especially, in a logarithmic time, chaos 
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becomes hidden [12]. Chaos, as an intrinsic nature of a given system, should 
not be affected by the choice of a coordinate gauge, and a chaos indicator, 
LE, should be defined as invariant under spacetime transformations. It means 
that the LE in general relativity should be expressed as a physical or so called 
proper quantity but not a coordinate quantity [16- 19]. 

Now we consider a particle, called "observer", moving along an orbit (not 
necessary to be a geodesic) in the spacetime with a metric ds 2 = g a pdx a dx 13 . 
In this letter Greek letters run from to 3 and Latin letters from 1 to 3. The 
observer can determine if his motion would be chaotic by observing whether 
the proper distances from his neighbors are increasing exponentially or not 
with his proper time. Both the distances and time are observables and should 
not depend on the choice of a coordinate system. Then we can apply the theory 
of observation in general relativity [16] to explore the dynamical behavior of 
the observer. 

Fig. 1 illustrates the trajectories of the observer and its neighboring particle 
called "neighbor" . The initial condition of the neighbor is arbitrarily chosen as 
long as its 4-dimensional distance from the observer is small enough to assure 
that the neighbor is approximately in the tangent space of the observer. In 
an arbitrary spacetime coordinate system x a , one can derive the equations 
of motion of the observer and its neighbor. Here we have to notice that the 
independent variable of the equations has to be the coordinate time t rather 
than their proper times r because the two particles have different proper 
times but one unique independent variable has to be adopted when integrating 
numerically their equations of motion together. At the coordinate time t the 
observer arrives at the point O with the coordinate x a and 4- velocity U a , and 
the corresponding proper time of the observer is r. At the same coordinate 
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time t, its neighbor reaches the point P with the coordinate y a along another 
orbit. A displacement vector 5x a = y a — x a from O to P should be projected 
into the local space of the observer. The space projection operator of the 
observer is constructed as h a ^ = g af3 + c~ 2 [/ a £/ /3 (c represents the velocity of 



light) [16]. OP' represents the projected vector <5x" — hgSx 13 , and its length 



\OP'\\ is 



AL(r) = ^g a(S 5xl5xi = ^h a(S 5x a 5x^ (5) 

Here g a p is calculated at x a . AL(r) is the proper distance to the neighbor 
observed by the observer at his time r and it is a scalar. Hence the maximum 
LE in general relativity is defined as 

1 , AL(t) 
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Let S r be the 3-dimensional subspace of the tangent space of the observer 
at O, which is orthogonal to the 4-velocity U a . It is the point P 1 but not P 
in S r as long as P' is close enough to the observer O. This tells us that Xr 
is coordinate invariant. The next is an argument for this point. Let us carry 
out a time transformation t — > rj. For convenience, we express the projection 
operation as OP' = HO?, where H is the space projection operator and 
= 0, where a subscript is added for if to describe the 4-velocity. Assume 
that the observer is located at the point O at r\ that corresponds to the old 
coordinate time t, then the neighbor at r\ would situate at the point Q but 
not necessary at P. Certainly we have to keep both P and Q inside an e 
neighborhood of O (see Fig.l) and e is considered as a very small quantity. In 
this case we have T? P = T? + 0(e). Furthermore, we have 0$ = UP — Qp 
and P$ = Arplfp + O(Arp) = keif P + 0(e 2 ), where k is a constant and 
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Fig. 1. The trajectories of the observer and its neighbor in spacetime. 

the proper time interval of the neighbor between P and Q, Arp, is in the 
magnitude of e. Hence, we get the relation HO($ = HOP + 0(e 2 ). This shows 
that \r is invariant with time transformations. 

As far as the numerical implementation of the computation of \r is concerned, 
the following notes are worth noticing. (1) The observer and the neighbor have 
different proper times, so a coordinate time t should be adopted as the inde- 
pendent variable. The equations of the motion of the particles should be trans- 
formed to use t as the independent variable. The variables to be computed 
step by step are the space coordinates and velocities (with respect to t) of the 
observer and its neighbor, and the proper time r and dr/dt of the observer. 
In total there are 14 variables. (2) As is well known, renormalization after a 
certain time interval is essential in this procedure to keep the distance be- 



6 



tween the observer and its neighbor small enough. The renormalization must 
be proceeded in the phase space though our LE is calculated in the configura- 
tion space. (3) One important difference during renormalization between the 
relativistic and classical cases must be noticed. Let E t be the local 3-space in 
which all the points have the same coordinate time t. It is evident that the 
renormalization should be done in H t otherwise the computation will commit 
an error. In Fig.l P is pulled back to M and Uti = (AL(0) / AL(t))oP . The 
next integration step for the neighbor will start from the point M. It is obvious 

that UP is located in T, t because both points O and P are in E t . On the other 
> 

hand OP is in the 3-space E T and a pull back from P' to S as renormalization 
is not correct. Our practice has proven this argument. 

In order to check the validity of our scheme for calculating the relativistic LE, 
we are going to reexamine the core-shell system studied by Vieira & Letelier[6]. 
In Schwarzschild coordinates (t,r,0,<f>), the 4-metric for this system is of the 
form 

ds 2 = -(1 - 2 -)e A dt 2 + e B - A [(l - ^dr 2 + r 2 d6 2 } 

+e _ V sin 2 6d<f) 2 , (7) 

where A and B are functions of r and 6 only. Here we adopt nondimensional 
variables and take c = 1. The metric does not explicitly depend on t and <fi 
and has an energy constant E and an angular momentum constant L, so test 
particles in free fall are actually in a system with two degrees of freedom with 
an integral U a U a = —1. When A = B = 0, Eq.(7) represents the Schwarzschild 
spacetime in which test particles in free fall move in regular orbits due to 
integrability of the system. For simplicity we discuss the dynamical behavior 
of geodesies in a black hole plus a dipolar shell, and let 



7 



A = 2<jfiv, 

B = l0 + Aav-a 2 [fi 2 (l-v 2 )+v% 



(8) 



where \i = r — \ and v = cos 9. The Newtonian limit of the model re- 
sembles the Stark problem[20,21], which is fully integrable. However, as far 
as the relativistic model is concerned, Vieira & Letelier[6] and Saa & Ve- 
negeroles[7] have demonstrated strong chaos in Weyl coordinates using the 
Poincare surface of section, and Gueron & Letelier[8] estimated its maximum 
LE Aiv = (3.2 ± 0.4) x KT 4 with the classical definition of LE (see Eq.(l)) in 
prolate spheroidal coordinates. 
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Fig. 2. The maximum Newtonian Lyapunov exponent, A at, with the time variable 
rj. It becomes one order of magnitude smaller after the time transformation (see 
Eq.(9)). 

We choose parameters as E = 0.975, L = 3.8, a = 2.5 x 10 -4 , and 7 = c 2 , 
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proper time x (10 5 ) 

Fig. 3. The maximum relativistic Lyapunov exponent, A_r, with the proper time r 
of the observer. It remains invariant after the time transformation (see Eq.(9)). r 
reaches about 9.168 x 10 5 when 77 runs 10 7 . 

and the initial conditions of the observer as r — 32, 9 — |, </> — 0, r = 0, 
and 9 from U a U a = —1. As to its neighbor, an initial separation Ar = — 1CT 8 
is adopted, regarded as the best choice [3], and the others remain the same 
as the observer's except 9. The values of E and L are carefully chosen to 
assure that the trajectories of the observer and its neighbors are bounded in a 
compact region. We integrate the geodesic equations using two integrators for 
comparison, Runge-Kutta -Fehlberg 7(8) and the 12th-order Adams-Cowell 
method, with a coordinate time step 0.01. When t reaches 10 6 , we find \n = 
\' N = (2.2 ± 0.2) x 10~ 4 , which is close to the result of [8] as its integration 
time amounts to 10 5 . Meanwhile, we get \ R = (2.8 ± 0.2) x 10" 4 by Eq.(6). 
One may notice that \n and Xr are in the same magnitude. This is because 
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the Euclidian distance d'(t) and the invariant proper distance AL differ not 
very much and so do the coordinate time t and the proper time r. In fact, r 
runs about 9.174 x 10 5 when t passes through 10 6 . A stronger gravitation by 
decreasing the angular momentum L will increase the difference between A at 
and X R , but the smallest stable circular orbit in Schwarzschild spacetime is 
r = 6. 

To verify the invariance of \r, we do a time transformation 

t -> n = lOt + r 2 /2. (9) 



We obtain = (2.6 ± 0.2) x 10 5 in the coordinate time 77, while \r retains 
the original value (see Fig. 2 and Fig. 3). 

As a further experiment, we slightly change the model to put B/2 = A = 
2o\iv. This system is still chaotic in Schwarzschild coordinates (t,r,9,(fi). To 
remove the singularity of this metric at the horizon, we go to Lemaitre coor- 
dinates (T,R,e,4>)[22] as 

T = /t(t + 2v / 27 + 21n| v ^~ ^ |), 

V r + v2 

^f^f + T, (10) 

where the positive constant k may be chosen arbitrarily. Our numerical tests 
display that the larger k is, the smaller A jv becomes. Particularly, the Lemaitre 
coordinates hide chaos for sufficient large k if \n is adopted as a chaotic index. 
However, our \ R does not vary with the parameter k. 

We would like to emphasize again that the computation of \r can be applied, 
whether the particles move along geodesies or not. The theory of observa- 
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tion in general relativity is not relevant to the 4-acceleration of the observer. 
For example, this method can be used to identify chaos in compact binary 
systems [23, 24]. 

In this letter we concentrate on calculating an invariant maximum LE in gen- 
eral relativity, but this discussion can be easily extended to find all the rel- 
ativistic LEs numerically with the technique proposed by Benettin et al.[25] 
(also see [26]). They suggest choosing an initial set of orthonormal tangent 
vectors as a base of the tangent space of the observer, then compute the evo- 
lution of the volume determined by these vectors. To avoid two vectors getting 
close to each other under evolution they use Gram-Schmidt procedure to or- 
thonormalize the base after each time step. The only change in a relativistic 
model is in computing the evolution of a tangent vector, which can be realized 
by the technique in this letter. As to the inner product in the ort honor maliz- 
ing procedure and the norm computation, one has to use the Riemanian inner 
product in place of the Euclidean one. 

We are very grateful to Dr. Xin-lian Luo of Nanjing University for his helpful 
discussion. This research is supported by the Natural Science Foundation of 
China under contract Nos. 10233020 and 10173007. 
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